Data Analysis: Out of the Lab and into the Streets?

Introduction

This document includes the statistical analyses reported in the paper Dablander, F.^⭑, Sachisthal, M.^⭑, & Aron, A. (under review) Out of the Labs and into the Streets? Effects of Climate Protests by Environmental Scientists. Royal Society Open Science.

Data preparation

We prepare the data and then test the above hypotheses.

library(brms)
library(dplyr)
library(knitr)
library(tidyr)
library(readr)
library(scales)
library(ggplot2)
library(corrplot)
library(gridExtra)
library(BayesFactor)
library(marginaleffects)
source('helpers.R')

df_text <- read_csv('data/dat_text.csv')
df <- read_csv('data/dat_num_all.csv')
demographics <- read.csv('data/demographics_all.csv')
df <- df %>%
  rename(Submission.id = SESSION_ID) %>% 
  left_join(
    demographics %>%
      rename(political_affiliation = U.s..political.affiliation) %>% 
      select(Submission.id, Age, Sex, political_affiliation),
    by = 'Submission.id'
  ) %>% 
  mutate(
    Protest = case_when(
      grepl('Legal', condition) ~ 'Legal March',
      grepl('CD', condition) ~ 'Civil Disobedience'
    ),
    Protest = factor(Protest, levels = c('Civil Disobedience', 'Legal March')),

    Engagement = case_when(
      grepl('None', condition) ~ 'Control',
      grepl('Endorse', condition) ~ 'Support',
      grepl('Join', condition) ~ 'Engaged'
    ),
    # Engagement = factor(Engagement, levels = c('Engaged', 'Support', 'Control')),
    Engagement = factor(Engagement, levels = c('Control', 'Support', 'Engaged')),

    PolicySupport = as.numeric(PolicySupport),
    ActivistSupport = as.numeric(ActivistSupport),
    ScienceCredibility = as.numeric(ScienceCredibility),
    SourceCredibility = as.numeric(SourceCredibility),
    Radical = as.numeric(Radical),
    Donation = as.numeric(Donation_1),

    PoliticalAffiliation = as.numeric(PoliticalAffiliation)
  ) %>% 
  filter(
    # !(political_affiliation %in% c('CONSENT_REVOKED', 'DATA_EXPIRED')),
    Progress == 100 # Remove 10 people who did not finish the survey
  ) %>% 
  mutate(
    political_affiliation = ifelse(
      PoliticalAffiliation %in% c(1, 2), 'Democrat',
      ifelse(
        PoliticalAffiliation %in% c(4, 5), 'Republican', 'Independent'
      )
    ),
    Politics = factor(political_affiliation, levels = c('Democrat', 'Republican', 'Independent'))
  )

contrasts(df$Protest) <- contr.sum(levels(df$Protest))
contrasts(df$Engagement) <- contr.sum(levels(df$Engagement))
contrasts(df$Politics) <- contr.sum(levels(df$Politics))

nrow(df)

## [1] 3149

We remove those participants that incorrectly answer the attention check question.

df$AC_Protest %>% table

## .
##    1    2    3 
## 1755 1350   44

df_text$AC_Protest %>% table

## .
##            Blockade with arrests Peaceful protest without arrests 
##                             1350                             1754 
##             There was no protest 
##                               44

df <- df %>% 
  filter(
    (Protest == 'Civil Disobedience' & AC_Protest == '2') |
    (Protest == 'Legal March'        & AC_Protest == '1')
  )

We are left with 2,875 responses.

nrow(df)

## [1] 2875

We remove the users that responded incoherently or inappropriately to the open text question about what the article was about. This leaves us with 2,856 responses.

df_open <- read.csv('data/openAnswers_coded.csv') %>% filter(Usable == 1)

# Additional check from second coder
df_check <- df_open %>% filter(Check == 1)
remove_responseids <- df_check$ResponseId[-c(4, 6)]
df <- df %>%
  left_join(df_open, by = 'ResponseId') %>% 
  filter(Usable == 1, !(ResponseId %in% remove_responseids))

nrow(df)

## [1] 2856

# Prepare sub-set of data for later analyses
df_march <- df %>% 
  filter(Protest == 'Legal March') %>% 
  mutate(Y = Radical)

df_civil <- df %>% 
  filter(Protest == 'Civil Disobedience') %>% 
  mutate(Y = Radical)

df_source <- df %>% 
  filter(!is.na(SourceCredibility)) %>% 
  mutate(
    Y = SourceCredibility,
    Engagement = factor(Engagement, levels = c('Support', 'Engaged'))
  )

contrasts(df_source$Engagement) <- contr.sum(levels(df_source$Engagement))

We fit all the models for the model-averaged exploratory analysis.

# Not included because it is ~50 MB large
fit_all <- readRDS('models/models_scaled.RDS')

# Not included because it is ~300 MB large
if (file.exists('models/fit_all_averaging_scaled.RDS')) {
  
  # For model-averaging later
  fit_all_averaging <- readRDS('models/fit_all_averaging_scaled.RDS')
  fit_policy <- fit_all_averaging$fit_policy
  fit_activist <- fit_all_averaging$fit_activist
  fit_radical <- fit_all_averaging$fit_radical
  fit_source <- fit_all_averaging$fit_source
  fit_science <- fit_all_averaging$fit_science
  
} else {
  
  set.seed(1)
  fit_policy <- get_logml_all_parallel(
    df %>% mutate(Y = PolicySupport), fit_all,
    iter = 4000, cores = 1, chains = 2
  )
  
  set.seed(1)
  fit_activist <- get_logml_all_parallel(
    df %>% mutate(Y = ActivistSupport), fit_all,
    iter = 4000, cores = 1, chains = 2
  )
  
  set.seed(1)
  fit_radical <- get_logml_all_parallel(
    df %>% mutate(Y = Radical), fit_all,
    iter = 4000, cores = 1, chains = 2
  )
  
  set.seed(1)
  fit_source <- get_logml_all_parallel(
    df_source, fit_all,
    iter = 4000, cores = 1, chains = 2, recompile = TRUE
  )
  
  set.seed(1)
  fit_science <- get_logml_all_parallel(
    df %>% mutate(Y = ScienceCredibility), fit_all,
    iter = 4000, cores = 1, chains = 2
  )
  
  fit_all_averaging <- list(
    'fit_policy' = fit_policy,
    'fit_activist' = fit_activist,
    'fit_radical' = fit_radical,
    'fit_source' = fit_source,
    'fit_science' = fit_science
  )
  
  saveRDS(fit_all_averaging, 'models/fit_all_averaging_scaled.RDS')
}

Confirmatory analysis

Policy support

fit <- fit_policy$fit_updated[[5]]
logmls_policy <- fit_policy$logmls

df_effects_policy <- get_all_predictions(
  fit, 'Policy support', variables = c('Engagement', 'Protest')
)

palette <- c("#E67C73", "#F2A489", "#F6E8A7", "#A3C4D4", "#5D92C7")
p1 <- create_plot(
  df, 'PolicySupport', 'Policy support',
  df_effects = df_effects_policy
)
p1

H1: Effect of protest on policy support

Moderate evidence against effect of protest.

h1_bfs <- get_bf_matrix(logmls_policy[c(1, 3)])
h1_bfs

##                H0      H1 H1r
## H0  vs  1.0000000 8.05561  NA
## H1  vs  0.1241371 1.00000  NA
## H1r vs         NA      NA  NA

H4: Effect of engagement on policy support

Strong evidence of no effect of engagement on policy support.

h4_bfs <- get_bf_matrix(
  c(logmls_policy[c(1, 2)],
    log(get_directed_bf(fit_policy$fit_updated[[2]])) + logmls_policy[2]
  )
)
h4_bfs

##                 H0       H1       H1r
## H0  vs  1.00000000 38.89893 43.953589
## H1  vs  0.02570765  1.00000  1.129944
## H1r vs  0.02275127  0.88500  1.000000

Activist support

fit <- fit_activist$fit_updated[[5]]
logmls_activist <- fit_activist$logmls

df_effects_activist <- get_all_predictions(
  fit, 'Activist support', variables = c('Engagement', 'Protest')
)

p2 <- create_plot(
  df, 'ActivistSupport',
  'Activist support',
  df_effects = df_effects_activist
)
p2

H2: Effect of protest on activist support

Overwhelming evidence in favour of a difference.

h2_bfs <- get_bf_matrix(
  c(logmls_activist[c(1, 3)],
    logmls_activist[3] + log(get_directed_bf(
      fit_activist$fit_updated[[3]], groups = 2,
      type = 'decreasing', pars = 'Protest'
    ))
  )
)
h2_bfs

##                   H0           H1          H1r
## H0  vs  1.000000e+00 1.745622e-15 8.728109e-16
## H1  vs  5.728617e+14 1.000000e+00 5.000000e-01
## H1r vs  1.145723e+15 2.000000e+00 1.000000e+00

H5: Effect of engagement on activist support

Strong evidence in favour of no effect, although it’s a bit lower compared to the directed hypothesis.

h5_bfs <- get_bf_matrix(
  c(logmls_activist[c(1, 2)],
    logmls_activist[2] + log(get_directed_bf(
      fit_activist$fit_updated[[2]]
    ))
  )
)
h5_bfs

##                 H0       H1      H1r
## H0  vs  1.00000000 14.43159 4.297033
## H1  vs  0.06929246  1.00000 0.297752
## H1r vs  0.23271871  3.35850 1.000000

Perceived radicalness of protest

fit <- fit_radical$fit_updated[[5]]
logmls_radical <- fit_radical$logmls

df_effects_radical <- get_all_predictions(
  fit, 'Perceived radicalness', variables = c('Engagement', 'Protest')
)

p3 <- create_plot(
  df, 'Radical',
  'Perceived radicalness',
  df_effects = df_effects_radical,
  lo = 'Not at all radical', hi = 'Extremely radical'
)
p3

H3: Effect of protest on perceived radicalness

Overwhelming effect of higher perceived radicalness for civil disobedience.

h3_bfs <- get_bf_matrix(
  c(logmls_radical[c(1, 3)],
    logmls_radical[3] + log(get_directed_bf(
      fit_radical$fit_updated[[3]], groups = 2,
      type = 'increasing', pars = 'Protest'
    ))
  )
)
h3_bfs

##                   H0           H1          H1r
## H0  vs  1.000000e+00 6.036804e-79 3.018402e-79
## H1  vs  1.656506e+78 1.000000e+00 5.000000e-01
## H1r vs  3.313011e+78 2.000000e+00 1.000000e+00

H6: Effect of engagement on perceived radicalness for legal march

Strong evidence evidence in favour of no effect of engagement.

h6_test <- get_logml_all(
  df_march,
  list('null' = fit_all[[1]], 'full' = fit_all[[2]])
)

## [1] 0.15
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5

# All marginal likelihoods
logmls_all <- c(
  h6_test$logmls,
  h6_test$logmls[2] + log(get_directed_bf(h6_test$fit_updated$full, type = 'decreasing'))
)

h6_bfs <- get_bf_matrix(logmls_all)
h6_bfs

##                 H0       H1      H1r
## H0  vs  1.00000000 4.301729 37.73446
## H1  vs  0.23246468 1.000000  8.77193
## H1r vs  0.02650097 0.114000  1.00000

H7: Effect of engagement on perceived radicalness for civil disobedience

Weak evidence evidence in favour of no effect of engagement.

h7_test <- get_logml_all(
  df_civil,
  list('null' = fit_all[[1]], 'full' = fit_all[[2]])
)

## [1] 0.15
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5
## Iteration: 6

# All marginal likelihoods
logmls_all <- c(
  h7_test$logmls,
  h7_test$logmls[2] + log(get_directed_bf(h7_test$fit_updated$full, type = 'decreasing'))
)

h7_bfs <- get_bf_matrix(logmls_all)
h7_bfs

##                H0        H1      H1r
## H0  vs  1.0000000 0.6879398 1.777622
## H1  vs  1.4536157 1.0000000 2.583979
## H1r vs  0.5625493 0.3870000 1.000000

Source credibility

We need to recompile the model because otherwise marginaleffects throws an error regarding the different factor levels for engagement being present.

# Is small enough to be included, but also quick enough to be estimated once
if (!file.exists('models/fit_source_full_final.RDS')) {
  fit <- brm(
    formula = 'Y ~ Engagement + Protest',
    data = df_source, family = brms::cumulative('probit'),
    prior = c(
      prior(student_t(1, 0, 0.15), class = b),
      prior(student_t(1, 0, 2.5), class = Intercept)
    ), cores = 4, chains = 4, iter = 1000, save_pars = save_pars(all = TRUE)
  )
  
  saveRDS(fit, 'models/fit_source_full_final.RDS')
  
} else {
  fit <- readRDS('models/fit_source_full_final.RDS')
}

df_effects_source <- get_all_predictions(
  fit, 'Source credibility', variables = c('Engagement', 'Protest')
) %>% 
  filter(Engagement != 'Control')

p4 <- create_plot(
  df_source, 'SourceCredibility',
  'Source credibility',
  df_effects = df_effects_source,
  lo = 'Not at all', hi = 'A great deal'
)
p4

H8: Effect of engagement on source credibility

Weak evidence in favour of no effect.

logmls_source <- fit_source$logmls
h8_bfs <- get_bf_matrix(logmls_source[c(1, 2)])
h8_bfs

##                H0       H1 H1r
## H0  vs  1.0000000 1.854321  NA
## H1  vs  0.5392809 1.000000  NA
## H1r vs         NA       NA  NA

Science credibility

fit <- fit_science$fit_updated[[5]]

df_effects_science <- get_all_predictions(
  fit, 'Science credibility', variables = c('Engagement', 'Protest')
)

p5 <- create_plot(
  df, 'ScienceCredibility',
  'Science credibility',
  df_effects = df_effects_science,
  lo = 'Not at all', hi = 'A great deal'
)
p5

H9: Effect of engagement on science credibility

Strong evidence for no effect of engagement.

logmls_science <- fit_science$logmls
h9_bfs <- get_bf_matrix(logmls_science[c(1, 2)])
h9_bfs

##                 H0      H1 H1r
## H0  vs  1.00000000 38.4238  NA
## H1  vs  0.02602553  1.0000  NA
## H1r vs          NA      NA  NA

Figure 1: Proportions and latent means

ggsave(
  'figures/Figure1.pdf',
  grid.arrange(
    p1, p2, p3, p5, p4,
    layout_matrix = rbind(
      c(1, 1, 2, 2),
      c(3, 3, 4, 4),
      c(NA, 5, 5, NA)
    )
  ),
  width = 11, height = 12
)

Table 1: Summary of confirmatory results

hyp_names <- c(
  'H1: No effect of protest form on policy support',
  'H2: Higher activist support for legal march',
  'H3: Higher perceived radicalness for civil disobedience',
  'H4: Higher policy support for scientist joining, then supporting, then for no engagement',
  'H5: Higher activist support for scientist joining, then supporting, then for no engagement',
  'H6: Legal march perceived as least radical when scientists join, then for supporting, then for no engagement',
  'H7: Civil disobedience perceived as least radical when scientists join, then for supporting, then for no engagement',
  'H8: No effect of type of engagement on source credibility',
  'H9: No effect of type of engagement on science credibility'
)
hyp_bfs <- c(
  h1_bfs[1, 2],
  h2_bfs[3, 1],
  h3_bfs[3, 1],
  h4_bfs[3, 1],
  h5_bfs[3, 1],
  h6_bfs[3, 1],
  h7_bfs[3, 1],
  h8_bfs[1, 2],
  h9_bfs[1, 2]
)

kable(
  cbind(hyp_names, round(hyp_bfs, 3)),
  col.names = c('Hypothesis', 'Bayes factor in favour of the hypothesis')
)

Hypothesis	Bayes factor in favour of the hypothesis
H1: No effect of protest form on policy support	8.056
H2: Higher activist support for legal march	1145723482010399
H3: Higher perceived radicalness for civil disobedience	3.31301116447911e+78
H4: Higher policy support for scientist joining, then supporting, then for no engagement	0.023
H5: Higher activist support for scientist joining, then supporting, then for no engagement	0.233
H6: Legal march perceived as least radical when scientists join, then for supporting, then for no engagement	0.027
H7: Civil disobedience perceived as least radical when scientists join, then for supporting, then for no engagement	0.563
H8: No effect of type of engagement on source credibility	1.854
H9: No effect of type of engagement on science credibility	38.424

Exploratory analysis

Here we show the results across political affiliation and report inclusion Bayes factors.

Policy support

df_effects_policy <- get_all_predictions(
  fit_policy$fit_updated[[19]], 'Policy support',
  variables = c('Engagement', 'Protest', 'Politics')
)

p7 <- create_plot(
  df, 'PolicySupport',
  'Policy support',
  'Expand offshore drilling',
  political = TRUE,
  df_effects = df_effects_policy
)
p7

ggsave('figures/FigureS4.pdf', p7, width = 9, height = 7)

inclusion_bfs_policy_support <- get_inclusion_bayes_factors(fit_policy)

kable(
  log(inclusion_bfs_policy_support),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	-3.3143482
Protest	-0.9983041
Politics	365.1042998
Engagement:Protest	-1.1661951
Engagement:Politics	-2.8801898
Protest:Politics	-2.2245190
Engagement:Protest:Politics	2.3180434

Activist support

df_effects_activist <- get_all_predictions(
  fit_activist$fit_updated[[19]], 'Activist support',
  variables = c('Engagement', 'Protest', 'Politics')
)

p8 <- create_plot(
  df, 'ActivistSupport', 'Activist support',
  political = TRUE, df_effects = df_effects_activist
)
p8

ggsave('figures/FigureS5.pdf', p8, width = 9, height = 7)

inclusion_bfs_activist_support <- get_inclusion_bayes_factors(fit_activist)

kable(
  log(inclusion_bfs_activist_support),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	-0.3266687
Protest	55.1913273
Politics	376.0767910
Engagement:Protest	0.7485845
Engagement:Politics	-3.6072906
Protest:Politics	-2.3995359
Engagement:Protest:Politics	-0.3930084

Perceived radicalness

df_effects_radical <- get_all_predictions(
  fit_radical$fit_updated[[19]], 'Perceived radicalness',
  variables = c('Engagement', 'Protest', 'Politics')
)

p9 <- create_plot(
  df, 'Radical', 'Perceived radicalness', 'of the protest',
  political = TRUE, df_effects = df_effects_radical,
  lo = 'Not at all radical', hi = 'Extremely radical'
)
p9

ggsave('figures/FigureS6.pdf', p9, width = 9, height = 7)

inclusion_bfs_radical <- get_inclusion_bayes_factors(fit_radical)

kable(
  log(inclusion_bfs_radical),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	1.701844
Protest	205.525590
Politics	108.389725
Engagement:Protest	-1.064140
Engagement:Politics	-3.230613
Protest:Politics	-2.129755
Engagement:Protest:Politics	-1.442445

Source credibility

# Is small enough to be included, but also quick enough to be estimated once
if (!file.exists('models/fit_source_politics.RDS')) {
  fit_source_full <- brm(
    formula = 'Y ~ Engagement + Protest + Politics',
    data = df_source, family = brms::cumulative('probit'),
    prior = c(
      prior(student_t(1, 0, 0.15), class = b),
      prior(student_t(3, 0, 2.5), class = Intercept)
    ), cores = 4, chains = 4, iter = 2000, save_pars = save_pars(all = TRUE)
  )
  
  saveRDS(fit_source_full, 'models/fit_source_politics.RDS')
  
} else {
  fit_source_full <- readRDS('models/fit_source_politics.RDS')
}

df_effects_source <- get_all_predictions(
  fit_source_full, 'Source credibility',
  variables = c('Engagement', 'Protest', 'Politics')
)

p10 <- create_plot(
  df %>% filter(Engagement != 'Control'),
  'SourceCredibility',
  'Source credibility',
  'Trust Dr. Fraser',
  political = TRUE, lo = 'Not at all', hi = 'A great deal',
  df_effects = df_effects_source
)
  
p10

ggsave('figures/FigureS7.pdf', p10, width = 9, height = 7)

inclusion_bfs_source <- get_inclusion_bayes_factors(fit_source)

kable(
  log(inclusion_bfs_source),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	0.2557725
Protest	-0.8414234
Politics	186.9059814
Engagement:Protest	-1.4305147
Engagement:Politics	-2.1577978
Protest:Politics	-1.0520638
Engagement:Protest:Politics	-0.6796638

Science credibility

df_effects_science <- get_all_predictions(
  fit_science$fit_updated[[19]], 'Science credibility',
  variables = c('Engagement', 'Protest', 'Politics')
)

p11 <- create_plot(
  df, 'ScienceCredibility',
  'Science credibility',
  'Trust environmental scientists',
  political = TRUE, df_effects = df_effects_science,
  lo = 'Not at all', hi = 'A great deal'
)
p11

ggsave('figures/FigureS8.pdf', p11, width = 9, height = 7)

inclusion_bfs_science <- get_inclusion_bayes_factors(fit_science)

kable(
  log(inclusion_bfs_science),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	-3.346037
Protest	-1.536171
Politics	421.843054
Engagement:Protest	-2.281769
Engagement:Politics	-3.599588
Protest:Politics	-1.760329
Engagement:Protest:Politics	-1.739139

Donation behaviour

fit_donation <- anovaBF(
  Donation ~ Protest + Engagement + Politics, data = df,
  rscaleEffects = c('Engagement' = 0.15, 'Protest' = 0.15, 'Politics' = 0.15)
)

if (!file.exists('models/df_effects_donation.csv')) {
  df_preds_donation <- get_df_preds_donation(fit_donation[18], df)
  df_effects_donation <- df_preds_donation %>% 
    group_by(Protest, Engagement, Politics) %>% 
    summarize(
      estimate = mean(ypred),
      conf.low = quantile(ypred, 0.025),
      conf.high = quantile(ypred, 0.975),
    )
  
  write.csv(df_effects_donation, 'models/df_effects_donation.csv', row.names = FALSE)
  
} else {
  df_effects_donation <- read.csv('models/df_effects_donation.csv')
}

df_effects_donation <- df_effects_donation %>% 
  mutate(
    Engagement = ifelse(Engagement == 'Support', 'Endorsed', Engagement),
    Engagement = factor(Engagement, labels = c('Control', 'Endorsed', 'Engaged'))
  )

cols <- c('#08519c', 'gray46', '#a50f15')
pdon <- ggplot(
  df_effects_donation, aes(x = Engagement, y = estimate, color = Politics)) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(width = 0.70), width = 0.50, size = 1
  ) +
  geom_point(size = 3, position = position_dodge(width = 0.70)) +
  theme_minimal() +
  labs(
    x = '',
    y = 'Predicted values',
    title = 'Donations across conditions'
  ) +
  scale_color_manual(values = cols) +
  scale_y_continuous(limits = c(10, 40)) +
  facet_wrap(~ Protest) +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 12),
    panel.grid = element_blank(),
    plot.title = element_text(size = 14, hjust = 0.50)
  ) +
  guides(fill = guide_legend(reverse = TRUE))

pdon

inclusion_bfs_donation <- c(
  # Main effects
  get_bf_inclusion_donation(fit_donation, 'Engagement')[1],
  get_bf_inclusion_donation(fit_donation, 'Protest')[1],
  get_bf_inclusion_donation(fit_donation, 'Politics')[1],
  
  # Two-way interaction effects
  get_bf_inclusion_donation(fit_donation, 'Engagement', varname2 = 'Protest')[2],
  get_bf_inclusion_donation(fit_donation, 'Engagement', varname2 = 'Politics')[2],
  get_bf_inclusion_donation(fit_donation, 'Protest', varname2 = 'Politics')[2],
  
  # Tree-way interaction effects
  get_bf_inclusion_donation(fit_donation, 'Engagement')[3]
)

names(inclusion_bfs_donation) <- names(inclusion_bfs_science)

kable(
  log(inclusion_bfs_donation),
  col.names = c('Effect', 'Log Inclusion Bayes factor')
)

Effect	Log Inclusion Bayes factor
Engagement	-2.390265
Protest	1.545690
Politics	29.764425
Engagement:Protest	-3.884956
Engagement:Politics	-5.542848
Protest:Politics	-4.276794
Engagement:Protest:Politics	-3.976229

Figure 2: Latent means across conditions

set.seed(1)
df_preds <- bind_rows(
  get_all_predictions(fit_policy$fit_updated[[19]], outcome = 'Policy support'),
  get_all_predictions(fit_activist$fit_updated[[19]], outcome = 'Activist support'),
  get_all_predictions(fit_radical$fit_updated[[19]], outcome = 'Perceived radicalness'),
  get_all_predictions(fit_science$fit_updated[[19]], outcome = 'Science credibility'),
  get_all_predictions(fit_source_full, outcome = 'Source credibility'),
  df_effects_donation %>% mutate(outcome = 'Donation')
) %>% 
  mutate(
    Politics = factor(Politics, levels = c('Democrat', 'Independent', 'Republican')),
    outcome = factor(
      outcome,
      levels = c(
        'Policy support', 'Activist support',
        'Perceived radicalness', 'Source credibility',
        'Science credibility', 'Donation'
      )
    ),
    panel = c(rep(1:4, each = 18), rep(5, 12), rep(16, 18))
  )

cols <- c('#08519c', 'gray46', '#a50f15')
ppost <- ggplot(
  df_preds, aes(x = Engagement, y = estimate, col = Politics, shape = Protest)
) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(width = 0.70), width = 0.50, size = 1
  ) +
  geom_point(size = 3, position = position_dodge(width = 0.70)) +
  theme_minimal() +
  labs(
    x = '',
    y = 'Latent mean',
    title = 'Latent means across outcomes and conditions'
  ) +
  scale_color_manual(values = cols) +
  facet_wrap(~ outcome, ncol = 2, scales = 'free_y') +
  theme(
    plot.title = element_text(size = 14, hjust = 0.50),
    strip.text = element_text(size = 10),
    legend.title = element_blank(),
    legend.position = 'top',
    legend.box = 'vertical',
    legend.spacing.y = unit(-0.25, 'cm')
  )

ggsave('figures/Figure2.pdf', ppost, width = 9, height = 9)
ppost

Table 2: Summary of exploratory results

varnames <- c(
  'Policy support',
  'Activist support',
  'Perceived radicalness',
  'Source credibility',
  'Science credibility',
  'Donation'
)

inclusion_bfs <- cbind(
  inclusion_bfs_policy_support,
  inclusion_bfs_activist_support,
  inclusion_bfs_radical,
  inclusion_bfs_source,
  inclusion_bfs_science,
  inclusion_bfs_donation
)

format_scientific <- function(x) {
  if (is.numeric(x) && x > 1000) {
    formatted <- formatC(x, format = "e", digits = 2)
    # Replace 'e+' with ' x 10^' for the desired format
    formatted <- gsub("e\\+0*", " x 10^", formatted)
    formatted <- gsub("e\\-", " x 10^-", formatted)
    return(formatted)
  }
  return(as.character(x))
}

colnames(inclusion_bfs) <- varnames
df_inclusion <- data.frame(round(inclusion_bfs, 3)) %>% tibble::rownames_to_column()
colnames(df_inclusion) <- c('Factor', str_replace_all(colnames(df_inclusion)[-1], '\\.', ' '))

# Apply the formatting function to the data frame
df_inclusion[, -1] <- as.data.frame(lapply(df_inclusion[, -1], function(column) {
  sapply(column, format_scientific)
}))

kable(df_inclusion)

Factor	Policy support	Activist support	Perceived radicalness	Source credibility	Science credibility	Donation
Engagement	0.036	0.721	5.484	1.291	0.035	0.092
Protest	0.369	9.32 x 10^23	1.81 x 10^89	0.431	0.215	4.691
Politics	3.65 x 10^158	2.13 x 10^163	1.18 x 10^47	1.49 x 10^81	1.60 x 10^183	8.44 x 10^12
Engagement:Protest	0.312	2.114	0.345	0.239	0.102	0.021
Engagement:Politics	0.056	0.027	0.04	0.116	0.027	0.004
Protest:Politics	0.108	0.091	0.119	0.349	0.172	0.014
Engagement:Protest:Politics	10.156	0.675	0.236	0.507	0.176	0.019

Supplementary Analysis

Sensitivity analysis: Confirmatory results

prior_widths <- seq(0.10, 0.50, 0.02) / 2

if (!file.exists('results/sensitivity_confirmatory.RDS')) {
    
  # Policy support
  h1_sens <- run_sensitivity(
    df %>% mutate(Y = PolicySupport),
    list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
    prior_widths
  )
  
  h4_sens <- run_sensitivity(
    df %>% mutate(Y = PolicySupport),
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths, includes_ordered = TRUE,
    groups = 3, pars = 'Engagement', type = 'increasing'
  )
  
  # Activist support
  h2_sens <- run_sensitivity(
    df %>% mutate(Y = ActivistSupport),
    list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
    prior_widths
  )
  
  h5_sens <- run_sensitivity(
    df %>% mutate(Y = ActivistSupport),
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths,
    includes_ordered = TRUE,
    groups = 3, pars = 'Engagement', type = 'increasing'
  )
  
  # Perceived radicalness
  h3_sens <- run_sensitivity(
    df %>% mutate(Y = Radical),
    list('null' = fit_all[[1]], 'full' = fit_all[[3]]),
    prior_widths
  )
  
  h6_sens <- run_sensitivity(
    df_march,
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths,
    groups = 3, pars = 'Engagement', type = 'decreasing'
  )
  
  h7_sens <- run_sensitivity(
    df_civil,
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths,
    groups = 3, pars = 'Engagement', type = 'decreasing'
  )
  
  # Source credibility
  h8_sens <- run_sensitivity(
    df_source,
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths
  )
  
  # Science credibility
  h9_sens <- run_sensitivity(
    df %>% mutate(Y = ScienceCredibility),
    list('null' = fit_all[[1]], 'full' = fit_all[[2]]),
    prior_widths
  )
  
  sens_all <- list(
    'h1_sens' = h1_sens,
    'h2_sens' = h2_sens,
    'h3_sens' = h3_sens,
    'h4_sens' = h4_sens,
    'h5_sens' = h5_sens,
    'h6_sens' = h6_sens,
    'h7_sens' = h7_sens,
    'h8_sens' = h8_sens,
    'h9_sens' = h9_sens
  )
  
  saveRDS(sens_all, 'results/sensitivity_confirmatory.RDS')
} else {
  sens_all <- readRDS('results/sensitivity_confirmatory.RDS')
}

Figure S2: Sensitivity analysis of confirmatory results

get_bf01 <- function(bfmat) bfmat[1, 2]
sapply(sens_all$h1_sens, get_bf01)

##  [1]  3.086449  3.548767  4.045935  4.481988  5.049997  5.524717  6.076725
##  [8]  6.485455  7.023044  7.590662  8.184941  8.699794  9.235297  9.596410
## [15] 10.182637 10.749493 11.251794 11.782758 12.399700 12.952234 13.457980

bfs <- do.call('c', lapply(names(sens_all), function(name) {
  sapply(sens_all[[name]], get_bf01)
}))

df_sens <- data.frame(
  type = rep(paste0('H[', seq(9), ']'), each = length(prior_widths)),
  prior_width = rep(prior_widths, 9),
  bfs = bfs
)

custom_label <- function(x) {
  ifelse(
    abs(x) > 500 | abs(x) < 0.002,
    formatC(x, format = 'e', digits = 1),
    format(x, scientific = FALSE)
  )
}

psens <- ggplot(df_sens, aes(x = prior_width * 2, y = bfs)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ type, scales = 'free_y', labeller = label_parsed) +
  labs(
    y = expression(BF[0*1]),
    x = 'Prior width',
    title = 'Bayes factor sensitivity analysis'
  ) +
  scale_y_continuous(labels = custom_label) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 12),
    plot.title = element_text(size = 14, hjust = 0.50)
  )

psens

ggsave('figures/FigureS2.pdf', psens, width = 9, height = 9)

Figure S3: Sensitivity analysis of exploratory results

prior_widths <- seq(0.10, 0.50, 0.04) / 2

policy_sens <- run_sensitivity_inclusion(
  df %>% mutate(Y = PolicySupport),
  fit_all,
  prior_widths,
  filename = 'results/policy_sens_incl.RDS'
)

activist_sens <- run_sensitivity_inclusion(
  df %>% mutate(Y = ActivistSupport),
  fit_all,
  prior_widths,
  filename = 'results/activist_sens_incl.RDS'
)

radical_sens <- run_sensitivity_inclusion(
  df %>% mutate(Y = Radical),
  fit_all,
  prior_widths,
  filename = 'results/radical_sens_incl.RDS'
)

science_sens <- run_sensitivity_inclusion(
  df %>% mutate(Y = ScienceCredibility),
  fit_all,
  prior_widths,
  filename = 'results/science_sens_incl.RDS'
)

source_sens <- run_sensitivity_inclusion(
  df %>% mutate(Y = SourceCredibility),
  fit_all,
  prior_widths,
  filename = 'results/source_sens_incl.RDS'
)

# For donation behaviour
donation_sens <- run_sensitivity_inclusion_donation(
  df,
  prior_widths,
  filename = 'results/source_donation_incl.RDS'
)

all_sens <- rbind(
  data.frame(policy_sens) %>% mutate(outcome = 'Policy support'),
  data.frame(activist_sens) %>% mutate(outcome = 'Activist support'),
  data.frame(radical_sens) %>% mutate(outcome = 'Perceived radicalness'),
  data.frame(science_sens) %>% mutate(outcome = 'Science credibility'),
  data.frame(source_sens) %>% mutate(outcome = 'Source credibility'),
  data.frame(donation_sens) %>% mutate(outcome = 'Donation')
)

df_sens_all <- all_sens %>%
  pivot_longer(cols = -c(prior_widths, outcome)) %>% 
  mutate(
    name = gsub('\\.', ' x ', name),
    name = factor(
      name,
      levels = c(
        'Engagement', 'Protest', 'Politics',
        'Engagement x Protest', 'Engagement x Politics',
        'Protest x Politics', 'Engagement x Protest x Politics'
      )
    ),
    outcome = factor(
      outcome,
      levels = c(
        'Policy support', 'Activist support',
        'Perceived radicalness', 'Source credibility',
        'Science credibility', 'Donation'
      )
    )
  )

cols <- c('#D55E00', '#A65628', '#029E73', '#56B4E9', '#0072B2', '#CC79A7')
psens2 <- ggplot(df_sens_all, aes(x = prior_widths * 2, y = log(value), color = outcome)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ name, scales = 'free_y') +
  labs(
    y = expression('Log Inclusion' ~ BF[1*0]),
    x = 'Prior width',
    title = 'Sensitivity analysis for exploratory inclusion Bayes factors'
  ) +
  # scale_y_continuous(labels = custom_label) +
  scale_color_manual(values = cols) +
  theme_minimal() +
  theme(
    legend.position = 'top',
    legend.title = element_blank(),
    strip.text = element_text(size = 12),
    plot.title = element_text(size = 14, hjust = 0.50),
    legend.direction = "horizontal",
    legend.box = "horizontal",
    legend.box.just = "center"
  ) + 
  guides(colour = guide_legend(nrow = 2))

psens2

ggsave('figures/FigureS3.pdf', psens2, width = 9, height = 9)

Supplementary Figure 9: Correlations between measures

df_cors <- df %>% 
  select(
    PolicySupport, ActivistSupport, SourceCredibility, ScienceCredibility,
    Radical, Donation, Gender, Age.x, Income, EducationLevel, PoliticalAffiliation
  ) %>% 
  mutate_all(as.numeric)

cormat <- cor(df_cors, method = 'kendall', use = 'pairwise.complete.obs')
diag(cormat) <- NA
rownames(cormat) <- colnames(cormat) <- c(
  'Policy support', 'Activist support', 'Source credibility', 'Science credibility',
  'Perceived radicalness', 'Donation', 'Gender', 'Age', 'Income', 'Educational level',
  'Political affiliation'
)

pdf('figures/FigureS9.pdf', width = 8, height = 8)
corrplot(
  cormat, method = 'color', number.cex = 0.80,
  addCoef.col = 'black', type = 'upper',
  tl.cex = 0.80, addrect = 20, tl.col = 'black',
  na.label = ' ',
  is.corr = TRUE,
  number.digits = 2, col.lim = c(-0.80, 0.80)
)
dev.off()

## quartz_off_screen 
##                 2

corrplot(
  cormat, method = 'color', number.cex = 0.80,
  addCoef.col = 'black', type = 'upper',
  tl.cex = 0.80, addrect = 20, tl.col = 'black',
  na.label = ' ',
  is.corr = TRUE,
  number.digits = 2, col.lim = c(-0.80, 0.80)
)

Session info

sessionInfo()

## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/Amsterdam
## tzcode source: internal
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] doParallel_1.0.17      iterators_1.0.14       foreach_1.5.2         
##  [4] stringr_1.5.1          marginaleffects_0.25.0 BayesFactor_0.9.12-4.7
##  [7] Matrix_1.6-5           coda_0.19-4.1          gridExtra_2.3         
## [10] corrplot_0.92          ggplot2_3.5.1          scales_1.3.0          
## [13] readr_2.1.5            tidyr_1.3.1            knitr_1.45            
## [16] dplyr_1.1.4            brms_2.21.0            Rcpp_1.0.12           
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1     farver_2.1.1         Rmpfr_0.9-5         
##  [4] loo_2.7.0            fastmap_1.1.1        tensorA_0.36.2.1    
##  [7] digest_0.6.35        estimability_1.5.1   lifecycle_1.0.4     
## [10] StanHeaders_2.32.6   magrittr_2.0.3       posterior_1.5.0     
## [13] compiler_4.3.3       rlang_1.1.3          sass_0.4.9          
## [16] tools_4.3.3          utf8_1.2.4           yaml_2.3.8          
## [19] collapse_2.1.0       data.table_1.17.0    labeling_0.4.3      
## [22] bridgesampling_1.1-2 bit_4.0.5            pkgbuild_1.4.4      
## [25] curl_5.2.1           abind_1.4-5          withr_3.0.0         
## [28] purrr_1.0.2          grid_4.3.3           stats4_4.3.3        
## [31] fansi_1.0.6          xtable_1.8-4         colorspace_2.1-0    
## [34] inline_0.3.19        emmeans_1.10.2       insight_1.1.0       
## [37] cli_3.6.2            mvtnorm_1.2-4        rmarkdown_2.27      
## [40] crayon_1.5.2         ragg_1.3.0           generics_0.1.3      
## [43] RcppParallel_5.1.7   rstudioapi_0.16.0    tzdb_0.4.0          
## [46] pbapply_1.7-2        cachem_1.0.8         rstan_2.32.6        
## [49] bayesplot_1.11.1     matrixStats_1.3.0    vctrs_0.6.5         
## [52] V8_4.4.2             jsonlite_1.8.8       hms_1.1.3           
## [55] bit64_4.0.5          systemfonts_1.0.6    jquerylib_0.1.4     
## [58] glue_1.7.0           codetools_0.2-19     distributional_0.4.0
## [61] stringi_1.8.3        gtable_0.3.5         QuickJSR_1.1.3      
## [64] gmp_0.7-4            munsell_0.5.1        tibble_3.2.1        
## [67] pillar_1.9.0         htmltools_0.5.8.1    Brobdingnag_1.2-9   
## [70] R6_2.5.1             textshaping_0.3.7    vroom_1.6.5         
## [73] evaluate_0.23        lattice_0.22-5       highr_0.10          
## [76] backports_1.4.1      bslib_0.7.0          rstantools_2.4.0    
## [79] MatrixModels_0.5-3   nlme_3.1-164         checkmate_2.3.1     
## [82] xfun_0.43            pkgconfig_2.0.3